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We study one-quasiproton excitations in the rare-earth region in the framework of the nuclear 
Density Functional Theory in the Skyrme-Hartree-Fock-Bogoliubov variant. The blocking prescrip- 
tion is implemented exactly, with the time-odd mean field fully taken into account. The equal filling 
approximation is compared with the exact blocking procedure. We show that both procedures 
are strictly equivalent when the time-odd channel is neglected, and discuss how nuclear alignment 
properties affect the time-odd fields. The impact of time-odd fields on calculated one-quasiproton 
bandhead energies is found to be rather small, of the order of 100-200 keV; hence, the equal filling 
approximation is sufficiently precise for most practical applications. The triaxial polarization of the 
core induced by the odd particle is studied. We also briefly discuss the occurrence of finite-size spin 
instabilities that are present in calculations for odd-mass nuclei when certain Skyrme functionals 
are employed. 



PACS numbers: 21.60.Jz, 21.10.Pc, 21.30.Fe, 27.70. +q 

I. INTRODUCTION 

The nuclear Density Functional Theory (dft) 0- 
3] plays a central role in the quest for a microscopic 
and quantitative description of atomic nuclei. The en- 
ergy functionals related to effective two-body density- 
dependent interactions are the main building blocks 
of the mean-field theory of the nucleus wherein the 
self-consistency is imposed through the Hartree-Fock- 
Bogoliubov (hfb) formalism. This framework has pro- 
vided a consistent description of a broad range of phe- 
nomena ranging from nuclear masses to collective exci- 
tations. Over the last few years, however, with the in- 
flux of high-quality experimental data on exotic nuclei, 
it has become evident that the standard local function- 
als (e.g., extended Skyrme functionals) are too restrictive 
when one is aiming at detailed quantitative description 
and extrapolability 0-0] • Consequently, various strate- 
gies have been devised to develop realistic nuclear energy 
density functionals (edf) @. These include: (i) the use 
of the density matrix expansion technique 0, llfj | to re- 
late the functional to low-momentum interactions; (ii) 
extending edf by adding higher-order terms in the local 
densities 11 1| : and (iii) improving spin and isospin prop- 
erties |12h15| . In any case, regardless of the strategy, the 
fine-tuning of the coupling constants of the functional to 
a suitably chosen set of experimental data is necessary to 
provide quality description (l6| . 

When aiming at spectroscopic-quality functionals @, 
the data coming from odd-mass nuclei are crucial: 
the energies, angular momenta, and parities of one- 
quasiparticle excitations provide us with basic knowledge 



about the underlying shell structure. Moreover, binding 
energies of odd- A systems are instrumental for determin- 
ing the magnitude of collective effects such as pairing. 
Theoretically, however, since nuclei with an odd number 
of particles have non-zero angular momentum ( J > 0), 
i.e., they are spin-polarized, their treatment is consider- 
ably more involved as compared to the J n —0 + ground- 
state (g.s.) configurations of doubly-even nuclei. 

Mathematically, the local edf is a time-even scalar 
constructed from various local densities and currents re- 
lated to particle and pairing density distributions [r7lll8j . 
The resulting mean field contains both time-even and 
time-odd terms. While the time-odd fields automati- 
cally vanish in the ground state of doubly-even nuclei, 
they are non-zero in J > configurations in which time- 
reversal symmetry is internally broken 0, Q . The time- 
odd fields have been investigated in the context of high- 
spin states [l9l - l22| , Gamow- Teller excitations [23[ , single- 
particle (s.p.) spectra d, [H, [HI , and collective dynamics 
[25T - l28| . The general consensus is that they can apprecia- 
bly impact the nuclear collective motion. On the other 
hand, our knowledge of the coupling constants charac- 
terizing individual time-odd fields is fairly limited, and 
the 

impact of those terms on nuclear ground-states still 
needs to be assessed. Conversely, one can ask whether 
experimental data on nuclear ground states can help con- 
strain the time-odd fields of the nuclear edf. 

There have been very few systematic theoretical stud- 
ies of one-quasiparticle states along isotopic or isotonic 
chains. Regional systematics of one-quasiparticle exci- 
tations, and their consequences on various observables 
in spherical and deformed nuclei, can be found in, e.g., 
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Refs. |29M33| (macroscopic-microscopic approach) and 
Refs. [jllMl34-l37| (nuclear DFT). The only global DFT 
study of ground state spin and parity for odd-mass nuclei 
is that of Bonneau et al. (38[. It is to be noted, however, 
that most of these studies were restricted in one way 
or another, e.g., by assuming axial symmetry, neglecting 
the time-odd fields, or doing an approximate treatment 
of blocking. The results of Refs. |g,LZl|38( clearly indicate 
that the currently used nuclear density functionals give a 
rather poor description of s.p. states, so it is imperative 
to evaluate the magnitude of the effects due to theoretical 
limitations and approximations. 

The goal of this study is to review the description of 
odd- mass nuclei in the framework of the nuclear DFT and 
assess the magnitude of time-odd polarizations through 
large-scale surveys. We compare various treatments of 
blocking, associated approximations, and resulting un- 
certainties. We discuss the choice of the orientation of 
the alignment vector, which is important for maintaining 
s.p. characteristics during the blocking procedure. We 
also assess the impact of the time-odd fields on binding 
energies of one-quasiparticle states and estimate the po- 
larization due to the axial symmetry breaking in certain 
orbits. 

This paper is organized as follows. Section |n] sum- 
marizes the main features of the nuclear Skyrme-DFT. 
We pay special attention to the treatment of odd-mass 
nuclei through the so-called blocking approximation and 
the Equal Filling Approximation (efa). In Sec. IIIII we 
present the details of the calculations and discuss vari- 
ous optimization techniques that enable large-scale cal- 
culations for odd-mass nuclei. The results are presented 
in Sec. lIVI We first compare the efa approximation with 
the exact blocking prescription. We estimate the effect 
of the time-odd fields on one-quasiparticle states in the 
rare-earth region and make selected comparisons with ex- 
periment. We also comment on the finite-size instabilities 
related to certain energy functionals that show up when 
studying polarized systems. Finally, the conclusions are 
contained in Sec. [Vj 



II. DFT TREATMENT OF 
ONE-QUASIPARTICLE STATES 



figuration representation, a basis of discrete states |n) is 
introduced, where n stands for all the s.p. quantum num- 
bers. The choice of one particular representation depends 
on the context. 

If |$) is the many-body state, the non-local density 
matrix in coordinate representation reads: 

p(m,rV) = ($10^,^1$), (1) 

where cJ. CT is a fermionic field operator creating a particle 
at position r with spin projection a and c ra is the corre- 
sponding annihilation operator. The field operators can 
be expressed in terms of the standard fermionic creation 
and annihilation operators and c n associated with the 
basis |n) gHH]: 

n 

Cra = E <p„(ra)c n . (2b) 

n 

Note that in this expression, </>„ {ra) and 0* (ra) are ma- 
trix elements of the basis transformation \ra) <H> \n): 
4>n(r<j) = (ra\n) and (/>*(rer) = (n\ra). They are there- 
fore complex s.p. wave functions dependent on the posi- 
tion vector r and spin coordinate a. The inverse relations 
are: 

4 = / <Pr^2,<f>n{ra)cl a , (3a) 
c n = d 3 r^2(/3* n (ra)c r(T . (3b) 

J a 

For complete bases, relations (|2all2bl) and (|3ati3bp allow us 
to express the relations between the two representations, 
p(ra,r'a') and p mn of the density matrix. 

The density matrix (fTJ) can be regarded as the ma- 
trix element of an operator p(ra, r'a') acting in the spin 
space. Any such operator can be expressed in terms of 
the Pauli matrices u and the identity matrix. This leads 
to a spin-scalar p(r,r') and a spin-vector field s(r,r'). 
These two fields are the fundamental building blocks of 
the nuclear dft. 



The nuclear dft in a Skyrme variant has been pre- 
sented in great detail in a number of articles [2|, |l8|, |39| . 
In the following we recall only the salient features of the 
theory that will be needed in this study. 



A. Representations of the Density Matrix 

The cornerstone of the nuclear dft is the general one- 
body density operator p. Two representations of the den- 
sity matrix are often considered. In the coordinate rep- 
resentation, the s.p. space is spanned by the continuous 
basis of states \ra) = \r) ® |<r) @, In tne con - 



B. Skyrme Energy Functional 

The contribution to the total energy of the system com- 
ing from the Skyrme interaction reads: 

^Skyrme = £ f j^ven) (j>) + ^(odd) (p) j ^ (4) 
t=0,l 

where t=0 and t—1 corresponds to isoscalar and isovec- 
tor components, respectively. In this paper, we do not 
consider proton- neutron mixing [l8| . Using the standard 
notation for the local densities and currents [13, [3 , the 
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part of the energy density that depends on time-even 
fields can be written as: 

ni cvcn \r) = c p tP 2 + cf? Pt Apt + ciPtn 

+ + C t VJ Pt V • Jti ( 5 ) 

while the part depending on the time-odd fields is: 

H t (odd) (r) = CI si + C t As s t ■ As t + Cjs t ■ T t 

+ CU\ + CT ] s t -(VAj t ). + C[s t -F t , (6) 

All densities and currents entering Eqs. ([5]) and ([5]) can 
be related to the particle density p{r, r'), spin density 
s(r,r'), and their derivatives [13, [ljj. In the present 
work, we do not consider tensor interactions and there- 
fore we set Cf = 0. 

Below, we discuss several versions of the functional, 
depending on how the time-odd coupling constants are 
determined: 

• Native version, which corresponds to all time-odd 
coupling constants being determined by the under- 
lying Skyrme interaction [l8[. 

• Gauge version, which corresponds to the subset 
of time-odd coupling constants being determined 
through the gauge- invariance conditions [H, l2lj |. 
namely, C\ = -Of, Cj = -C/, and CV = C? J , 
and all other time-odd coupling constants set to 
zero. 

• Landau version, which is based on the gauge version 
where the subset of time-odd coupling constants C t s 
and Cj are reset through the Landau parameters 
23]: 

.go = Ao (2C S + 2Cj(3p 2 Q /3 ) , 
9l = -2N Cjpp 2 /3 , 
g' = N a (2Ci+2C{f3p 2 /3 ), 
g[ = -2N CT(3p 2 /3 , 

where /3 = (3tt 2 /2) 2 / 3 , 1/Nq = TT 2 h 2 /2m*k F and, 
additionally, C t As = for t = 0, 1. Since the Lan- 
dau prescription only sets C t T , the gauge condition 
is broken since Cj ^ — C/ anymore. 

• Time-even version, in which all time-odd coupling 
constants in Eq. ([5]) are set equal to zero. 

C. HFB Method 

In the HFB theory, pairing correlations enter through 
the pairing tensor k defined in coordinate representation 
as: 

K (ra,r'a') = {$\c r , (7 ,c r(T \$), (8) 



(From a practical point of view, it is sometimes more 
advantageous to use the pairing density p |4ll |42T|.) 

The starting point of the hfb theory is to assume that 
the ground-state of an even-even nucleus is a vacuum for 
quasiparticle operators (/3„,/3j). The latter are obtained 
from single-particle operators (c„, c*) associated with the 
single-particle basis states |n) by the Bogoliubov trans- 
formation: 

P V = Y, U nu^ + V* v cl (9a) 

n 

Pt=J2 V ™ Cn + Un » C n- ( 9b ) 

71 

The matrices U and V are obtained from the HFB equa- 
tions: 

(ill -AJ( U v)-*( D v)- « 

where A is the chemical potential, h is the Hartree-Fock 
(hf) potential and A the pairing potential. (From a prac- 
tical point of view, it is sometimes more advantageous 

to use the pairing potential h [U H^].) The form of 
the hfb equations in coordinate space can be found in 
Refs. 

The density matrix and pairing tensor can be written 
as: 

Pmn={V*V T ) mn , (11a) 
n mn = (V*U T ) mn . (lib) 

The coordinate representation of the Bogoliubov trans- 
formation, 

ft,= f d 3 r5^{l7M*(r £ r)c r<7 + vW*(r(7)4 ff },(12a) 
Pt = f d 3 rY,{v (v) (r<T)c r<T + U^{ra)ci a ) ,(12b) 

can be expressed through lower and upper components 
of the quasi-particle (q.p.) wave functions: 

V^(ra) = J2K(r<j)V n „ (13a) 

71 

U^(r<T) = Y,^( r<T ) U ™- ( 13b ) 

n 

Finally, the density matrix and pairing tensor in coordi- 
nate space are: 

p{ra,r'a')= ^ V M *(r(x)V M (r'a'), (14a) 

0<B M <B max 

K{ra,r'a')= ^ V^'HE/WfrV). (14b) 

0<E„<E max 

It is assumed here that the q.p. continuum with E > —A 
has been discretized and all q.p. states with energy lower 
than some cut-off energy E cut are retained (see discussion 
in Ref. jij). 



4 



D. The Blocking Prescription and the Equal Filling 
Approximation 

In the hfb theory, the ground-state of an odd nucleus 
is a one quasiparticle excitation f}> with respect to the 
q.p. vacuum. In the configuration representation, the 
corresponding density matrix and pairing tensor are [45I — 
El: 

P%'%° = (V*V T ) mn + U m , U; w> - KVo^mo, (15a) 
«™r = (V*U T ) mn + U mfl0 V^ ~ V^ U nfi0 . (15b) 

In practice, one must adopt a prescription to be able to 
determine, at each iteration, the index /iq of the quasi- 
particle state to be blocked [49]. In the present study, 
this has been clone according to the recipe described in 
Ref. In the first step, the mean field Hamiltonian h 
is diagonalized: 

Jvp n = e n ip n . (16) 

Since in this work parity and y-signature are assumed 
to be self-consistent symmetries, every s.p. level e„ is 
uniquely identified by its position in a given parity and 
y-signature block. This unique identification allows to 
pin down the configuration of the blocking candidate no. 
To connect the s.p. state cp„ with a quasiparticle state 
to be blocked, we calculate at each iteration the overlap 
between tp no and both the upper component and the 
time-reversed lower component Vp of quasiparticle states 
around the Fermi level [5(|. The largest overlap in this 
set defines the index of the quasiparticle state to be 
blocked. In the beginning of the iterative process, s.p. 
states of a neighboring even-even nucleus are taken. 

Within the efa, the states /io and its time-reversal 
partner £lq enter the density matrix and pairing tensor 
with the same weights, which ensures time-reversal sym- 
metry and thereby degeneracy of /j,q and £lq [HI : 




(17b) 

The hfb equations are then solved by replacing (p B ' w , 
k B '^°) with (p EFA ^o ; k efa :W) ). p or the justification of the 
EFA ansatz by means of statistical density operators and 
for detailed discussion of the procedure involved, we refer 
the reader to Refs. [H,[5l|. 

In this work we point to another possible justification 
of the efa. We first notice that the time-even parts of the 
blocked density matrices given by Eqs. (|15a[) and (|15b[) 
are identical to the time-even parts of the density matri- 
ces in efa, Eqs. (|17al) and ()17b[) . Therefore, all time-even 



densities in Eq. ([5]) are exactly the same in both variants. 
Consequently, in the blocking and efa approximations, 
the time-even part of the functional (Sec. Ill B|) yields ex- 
actly the same self-consistent solution. This allows us to 
reinterpret efa density matrices as those corresponding 
to the time-even functional in which the time-odd po- 
larizations exerted by the odd particle are dynamically 
switched off. Of course, the blocking prescription and 
efa give exactly the same average values of all time-even 
observables (e.g., radii and multipole moments) but dif- 
fer in the average values of time-odd observables (e.g., 
spin alignments and magnetic moments). 

E. Blocking, Alignments, and Symmetries 

Although for the functionals restricted to time-even 
fields (or within efa) the time-reversed q.p states \fi) = 
/^|0) and |/2) = /3p|0) are exactly degenerate, this is not 
true any more in the general case. Here, the blocking 
prescription does depend on which of these two states, or 
which linear combination thereof, is used in Eqs. (|15a[) 
and (|15b[) . In order to discuss this point, we introduce 
here the notion of an "alispin" , which pertains to the 
unitary mixing of states and \p). This is in complete 
analogy with the standard notion of the isospin, which 
involves the unitary mixing of proton and neutron states. 

An alivector is defined as a set of two complex 
numbers a and b (\a\ 2 + \b\ 2 = 1): 

V W = ( I ) , (18) 

which corresponds to the linear combination of states \fi) 
and |/2): \v^) = a\fi) +b\p). Alivectors reside in SU(2) 
space; therefore the alirotation by an angle <p^ p ' is defined 
as: 

V G0' = e ^W») fay (1Q) 

where the alivectors of Pauli matrices are denoted by 
er^\ and o denotes the scalar product of alivectors. To 
recall that the alirotation pertains to a single pair of 
states, we use superscripts (/x) throughout. 

The blocked density matrix Eq. (|15a[) corresponding to 
the state reads: 

P*i a ' b) = Pmn ~ {\a\ 2 V nfl V^ + \b\ 2 V n ^ fi 

+a*bV n p,V£ lll + ab*V nfl V r * nfl } 
+ {\a\ 2 U^U m ^ + \b\ 2 U^U mfi 

+a*bU^U m ^ + ab*U r ^U mfi }. (20) 

If time-reversal symmetry is conserved, the different 
blocks of the Bogoliubov matrices are related: 

~V^ = V nn and V^ = V m , (21) 
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and pmn) = Pfh!n* ■ These relations lead to: 

n ( ab ) = n _ v v* +77* II (11) 

Hmn — Fran V n/J, ' mn > u nu u ran, K^^l 

or, equivalently 

Pmn — Pmn ~ VnpYmfi ^n^Pmfi- (23) 

Therefore, in this limit, the exact blocking density ma- 
trix becomes independent of the coefficients (a, b) of the 
mixing, i.e., it is an aliscalar. Since 



EFA 
Pmn 



0) 



. p (0,l) 
rrnn 



(24) 



the EFA density matrix also coincides with the exact 
blocking density matrix; hence, it is an aliscalar as well. 

In the general case where time-reversal symmetry is 
not dynamically conserved, however, the blocking den- 
sity matrix is not aliscalar and the energy of the sys- 
tem may change as a function of the mixing coefficients 
(a, b). To analyze what are the consequences of blocking 
different alirotated states V^o) , we introduce the (real) 
alignment vector — (p,\J\p) = — {p\J\p,) and the 

(complex) decoupling vector D^' — (fi\J\p.) = (fl\J\p)* 
52]. Together, they form the matrix elements of the 
alignment vector- alivector 



D 



(25) 



Expanding this operator (acting on SU(2) alistates) in 
the basis of Pauli matrices, we find: 



ft 
ft 
ft 



(26a) 
(26b) 
(26c) 



where indices k — 1,2,3 enumerate the components of 
alivectors. From these considerations it follows that the 
alignment vector-alivector J^ 1 , which corresponds to the 
alirotated pair (|19p . is obtained by: 



ft 
ft' 
ft' 



— exp(i(; 



'a ° S n ) 



ft 
ft 
ft 



(27) 



where are the standard spin-1 matrices 



, which 

are generators of rotation in the vector representation. 
This shows that the concept of alirotation (equivalent to 
changing the mixing of the blocked state) translates into 
a change in the alignment of the system. 

To illustrate how this works, let us examine the spe- 
cial case where the states \p) = \p y ) and \p) = \p, y ) are 
eigenstates of the R y signature operator. Since: 



R^Rj — ^ ^ £jjkRk ; 



(28) 



we can express the states 
eigenstates: 



P-y, 



and 



Py, 



in terms of R z 



\Py) 



\Py) 



1 
1 



V V2 



1 

1 



V2 



(29) 

which corresponds to the rotation of the system by the 
Euler angles (a, /3,7) = (0, n/2, ir/2). Consequently, the 
alivector V<>m) is the vector alirotated by the an- 

gles = (a, 0,7). Since (i) the blocked density ma- 
trix is not an aliscalar, and (ii) alirotations are induced 
by rotations of the coordinate system or, equivalently, 
a change of the symmetry operators used to label s.p. 
and q.p. states, we must conclude that the blocked den- 
sity matrix may depend on the choice of the symmetry 
operators that commute with the Hamiltonian 1 . 

More generally, since all alignment properties of the 
system are embedded in the vector-alivector operator J ', 
we also see that the alirotation of states \p) and \ corre- 
sponds to hfb states having different alignment vectors. 
Therefore, the latter can be used to tag blocked states. 
This is a very convenient method, which can be applied 
not only for the time-even version of the functional when 
the quasiparticle states \p) and \p,} are degenerate, but 
also in the case of internally broken time-reversal sym- 
metry. The key to our considerations of blocked states is 
the realization that blocking must depend on the orien- 
tation of the alignment vector with respect to the prin- 
cipal axes of the mass distribution. Therefore, the only 
rigorous way to proceed would be for each quasiparticle 
excitation to vary the orientation of the alignment vector 
with respect to the principal axes of the system, and re- 
tain the solution with the lowest energy [54[ . We give in 
Sec. IIV B 31 a pedagogical illustration of such anisotropy 
of blocking. 

In many practical applications, however, one chooses 
a fixed direction of alignment dictated by practical con- 
siderations. In particular, the identification of blocked 
single-particle states n and quasiparticle states /io is most 
conveniently carried out through the set of conserved 
quantum numbers characteristic of the problem. In all 
calculations performed in this work, nuclei are either axi- 
ally deformed or nearly axial, and they conserve reflection 
symmetry. The corresponding symmetry group is D^P; 
hence, signature r = ±i and parity n = ±1 are good 
quantum numbers. In hfodd signature is defined with 



1 Note that, in the particular case where the alivector is built from 
the eigenstates \fi v ) and \p, y ) of R y , the alirotation by (0, tt,0) is 
equivalent to the R y symmetry, and therefore leaves the sys- 
tem invariant. This operation corresponds to (\fJ. y ), \p. y )) — > 
~\t 1 y))- Therefore, in this particular case, blocking state 
\fj, y ) or state \p, y ) gives exactly the same energy, even though 
time-reversal symmetry is internally broken and the q.p. spectra 
do not exhibit the Kramers degeneracy. 
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respect to the y-axis of the reference frame [55|. In this 
way, the alignment vector is restricted to having only the 
y-component. To realize the three possible alignments of 
the angular momentum along the principal axes, it is suf- 
ficient to orient the longest, shortest, or intermediate axis 
along the y-axis. Since in most cases, the configurations 
analyzed in this study are axial, only two orientations 
suffice. We show in Sec. IIII Al below how to implement 
such a scenario. 

Equivalently, one could work with a good z-simplex 
basis such as in Ref. (56J. In that case, the default align- 
ment is along the z-axis, but results still depend on the 
orientation of the body. It is only by allowing the align- 
ment vector to cover the full solid angle that physical 
properties of the system would not depend on the choice 
of the basis used to describe the odd nucleus. 



III. METHOD OF CALCULATION AND 
OPTIMIZATION TECHNIQUES 

This section briefly describes the DFT solvers used in 
this work and discusses the choice of parameters entering 
our calculations. We also outline various optimization 
techniques that have been implemented by us to carry out 
large-scale dft calculations for one-quasiparticle states 
on leadership class computers. 



A. Numerical Parameters 

All calculations in this work are performed with the 
dft solvers hfbtho [13] and HFODD [5(| HE ■ Both 
codes solve the Skyrme hfb problem in the configuration 
space by means of the Harmonic Oscillator (ho) expan- 
sion technique. In hfbtho, the cylindrical HO basis is 
used, and both axial and time-reversal symmetries are 
imposed. This implies that the efa must be used for 
blocking calculations. The 3D solver hfodd employs the 
Cartesian HO basis and is symmetry-unrestricted. This 
unique feature of hfodd makes it a tool of choice for 
our study, since in the polarized nuclear configurations 
many self-consistent symmetries are usually broken. The 
blocking prescription is implemented exactly in hfodd 
with all the time-odd fields taken into account. Both 
codes have been benchmarked against one another and 
they yield the same results within a few eV for spherical 
or axially deformed even-even nuclei (59|. 

As already mentioned, all nuclei considered in this 
work are either axial or slightly triaxial, as well as 
reflection-symmetric. Therefore y-signature and parity 
are conserved and used in hfodd to tag q.p. and single- 
particle states. However, this implies that the total align- 
ment is confined to the y-axis. Since the latter is not 
the quantization axis, one can not easily associate the 
single-particle spin Q with the expectation value of the 
angular momentum: the situation is analogous to the col- 
lective rotation in high-spin physics. For the sake of iden- 



tification of deformed Nilsson orbitals, it is convenient, 
however, to reintroduce Q as a (nearly) good quantum 
number by orienting the angular momentum along the 
z-axis; the resulting alignment properties correspond to 
the limit of non-collective rotation. 

To this end, we need to associate the quantization axis 
with the symmetry axis of the nucleus. This can be 
achieved via a Euler rotation of the body-fixed frame 
(by a = it/2, f3 = n/2, 7 = 0) or by imposing con- 
straints on the expectation values of the quadrupole ten- 
sor Q2n- After testing these two options, we choose the 
Euler rotation: calculations for even-even nuclei are per- 
formed in the standard y-signature mode, then solutions 
are Euler-rotated and used to warm-start calculations for 
odd-nuclei. In this way, the Nilsson quantum number f2 
is computed from the expectation value of j y . This tech- 
nique turns out to be both stable and fast. Note that 
the energies of a given blocked state in the Euler-rotated 
case and original orientation are different, as discussed in 
Sec. Ill El Only a complete survey of all possible orien- 
tations of the alignment vector, which would be a major 
computational endeavour, could pin down the correct ori- 
entation. 

As is well known, calculations for deformed nuclei con- 
verge faster if the eigenstates are expanded on a stretched 
basis that follows the geometry of the nuclear density. 
Unfortunately, the stretched basis is not compatible with 
the Euler-rotation of the nucleus in space in HFODD. For 
that reason, all calculations presented in this work have 
been carried out in a full spherical basis of N osc =14 os- 
cillator shells (the number of basis states is N s = 680). 
This choice guarantees stability of results for the rel- 
atively modest deformations considered in this study. 
The oscillator frequency was fixed at 1.2 x Hujo [13] for 
huj = 4l/A 1 / 3 MeV. 

In this work, we use three commonly used Skyrme 
parametrizations: SIII SkP [U, and SLy4 (62[|. In 
the pairing channel, we employ the density-dependent 
delta interaction in the mixed variant 1631: 



V(r,r') = V 



Mr) 
2 Po 



(30) 



where po = 0.16fm _1 and Vo is the pairing strength (iden- 
tical for protons and neutrons) . Note that the use of such 
a zero-range interaction requires us to introduce a renor- 
malization (or regularization) procedure to avoid non- 
physical divergences [H], We employ the standard 
value i; ctlt =60MeV. 

For each Skyrme edf, the pairing strength Vo has been 
adjusted to reproduce the experimental proton odd-even 
mass difference in the deformed nucleus 162 Dy, Ap 3 " 1 = 
0.60 MeV. This choice has been motivated by the find- 
ings of Ref. [66| that by adjusting Vq to experimental 
data for a spherical semi-magic nucleus one underesti- 
mates pairing correlations in deformed systems having 
lower single-particle level density around the Fermi sur- 
face. Moreover, by considering the proton pairing gap, 
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one effectively takes into account the Coulomb contri- 
bution to pairing [67]]. The pairing strengths used in 
this work are V = -314.406 MeV, -297.303 MeV, and 
-2 49.059 Me V for SLy4, SIII, and SkP, respectively. In 
Sec uTEl we also consider SkO, V = -269.226 MeV, and 
SkM*, Vq = —297.875 MeV. 



B. Parallelization and Optimization 

The advent of Teraflop and Petaflop supercomputers 
enables large-scale surveys with symmetry-unrestricted 
dft solvers. To optimize resources, however, optimiza- 
tion of the production codes is required. Starting from 
the original published versions of hfodd and HFBTHO, 
we made a number of improvements. First of all, a par- 
allel interface using the standard Message Passing Inter- 
face (MPI) has been constructed to allow the automated 
distribution of calculations over several computing cores. 
Let us note that the standard nuclear dft calculations 
are "embarrassingly parallel". Indeed, solving the HFB 
equations for one nuclear configuration usually does not 
exceed a few hours on a standard desktop computer. 
Therefore, each computing core of a massively parallel 
system can process a single HFB task corresponding to 
a particular nucleonic configuration. Only in the limit 
of very large HO bases, or for dft solvers constructed 
in the coordinate space, does the parallelization of the 
solvers become necessary. The advantage of using mas- 
sively parallel architectures is that simultaneous calcula- 
tions of hundreds or thousands of different many-body 
configurations are possible in a very reasonable time. 
Such a strategy makes it possible to extract systematic 
trends, use standard statistical analysis toolboxes, and 
ultimately develop nuclear edfs of spectroscopic quality. 

The scaling of a dft problem with the number of pro- 
cessors also implies that a simple master-slave parallel ar- 
chitecture is sufficient for most applications, and this so- 
lution is adopted here. All calculations in this study were 
performed on the Cray XT3/XT4 Jaguar supercomputer 
at the National Center for Computational Science at the 
Oak Ridge National Laboratory and on the Cray XT3 
Franklin supercomputer at the National Energy Research 
Scientific Computing Center at the Lawrence Berkeley 
National Laboratory. Typical production runs have in- 
volved from 8,000 to 12,000 computing cores per run, and 
the typical calculation time was about 2 hours. 

The HFB equations represent a coupled system of 
non-linear equations for nucleonic densities. The non- 
linearity enters through the dependence of the mean 
fields on densities (self-consistency) . In order to diminish 
the number of iterations required to attain a given pre- 
cision, we employ the modified Broyden method [68, 69]. 
The method is based on the observation that the con- 
vergence of the HFB process stops when the characteris- 
tic variables in the problem, for example the density p, 
do not change any more from one iteration to the next: 
_ p(™)|| ~ o. In other words, the HFB equations 



can be viewed as a fixed-point problem, and iterations 
can therefore be optimized by employing a quasi-Newton 
method. It was shown that the computational cost (in 
units of number of iterations) could be reduced substan- 
tially, by a factor of 3 to 4. Our particular implemen- 
tation of the modified Broyden method was described in 
Ref. [Hi]. 

As mentioned earlier, hfodd solves the hfb problem 
by expanding eigenvectors in the Cartesian HO basis. In 
terms of speed, one of the major bottlenecks in hfodd 
is the diagonalization of the hfb matrix. The latter is 
carried out with the subroutine ZHEEVR of the LAPACK li- 
brary. We found that a significant gain in terms of speed, 
up to 30-40% for large HO bases with, e.g., iV osc =20 
shells, could be obtained by using machine-specific im- 
plementations of the BLAS and LAPACK libraries, such as 

ATLAS. 



C. Massively Parallel Calculations: Convergence 
Improvements 

Based on the experience gained in self-consistent cal- 
culations for odd- A nuclei, it appears that calculations 
involving blocking are always less stable than those per- 
formed for even-even nuclei. Apart from the specific issue 
related to finite-size instabilities addressed in Sec. IIV El 
below, these numerical instabilities are related to the 
need to select, at each iteration, the blocked q.p. state. 
The blocking procedure is outlined in Sec. Ill Dl and the 
detailed justification can be found in, e.g., Ref. The 
selection method implies that the blocked q.p. state may 
change from one iteration to the next, in particular at 
the beginning of the calculation. This numerical noise is 
the price to pay for the full self-consistency, and it ex- 
plains why small differences in the initial conditions can 
actually affect the convergence process. 

When only a few nuclei are considered, and a small 
number of blocked configurations near the ground-state 
is calculated, one can often find ways to converge calcu- 
lations, such as: (i) changing the linear-mixing param- 
eter of the self-consistent scheme; (ii) starting from the 
unblocked fully-paired state correspon ding to an odd av- 
erage particle number (false vacuum) [70j; (iii) starting 
from the even-even nucleus with one more particle for 
a particle-like blocked state and with one particle less 
for a hole-like blocked state, as implemented in Ref. [49[ 
(we also used this method in our calculation); (iv) us- 
ing different values of the linear-mixing parameter for 
time-even, time-odd, and/or pairing fields. Whenever a 
blocking calculation fails to converge, one may repeat it 
by using one or several of these tricks until a converged 
result is obtained. This is what was done in previous 
studies involving self-consistent calculations, and it was 
possible, because these studies were focused on ground- 
state properties and only a minimum number of different 
configurations was considered. 

In our case, however, we consider thousands of con- 
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TABLE I: Comparison of EFA (hfbtho) with exact blocking (hfodd) for four one-quasineutron configurations in 121 Sn. The 
time-odd fields are switched off. Shown are: the quasiparticle energy E qp , neutron chemical potential An, neutron pairing energy 
Ep air , average neutron pairing gap A n = Tr(Ap)/iV , total r.m.s. radius, axial quadrupole deformation /3, total quadrupole 
moment Qtot, kinetic energy E^in (for protons and neutrons), total spin-orbit energy Eso, direct Coulomb energy Edir, and 
total energy E to t- The last two lines show the hfodd alignments of the blocked quasiparticles: Jy was calculated in the 
non-collective orientation and J± in the collective orientation, see Discussion in Sec. IIV B 31 In the EFA total alignments are 
equal to zero by construction. The orbits are labelled by the 1^(0,^) quantum numbers. The SLy4 functional is used in the 
particle-hole channel and the density-dependent delta interaction with Vo= _ 285.634 MeV. The differences between hfbtho and 
hfodd results are highlighted in boldface. 



Exact State 


dz/2 


(l/2 + ) 


^3/2 


(3/2+) 




(5/2-) 


37/2(7/2+) 


Blocking 


EFA 


Exact 


EFA 


Exact 


EFA 


Exact 


EFA 


Exact 


E qp (MeV) 


1.0076 


1.0080 


1.1822 


1.1820 


1.4570 


1.4570 


2.2879 


2.2880 


Xn (MeV) 


-7.7496 


-7.7496 


-7.7288 


-7.7288 


-7.9834 


-7.9836 


-7.6371 


-7.6371 


E% air (MeV) 


-9.2949 


-9.2948 


-9.4411 


-9.4410 


-8.7145 


-8.7147 


-10.4041 


-10.4036 


A n (MeV) 


1.0575 


1.0575 


1.0667 


1.0667 


1.0395 


1.0395 


1.1206 


1.1206 


r.m.s (fm) 


4.6895 


4.6895 


4.6889 


4.6889 


4.6894 


4.6894 


4.6884 


4.6884 


P 


-0.0257 


-0.0257 


0.0131 


0.0131 


0.0099 


0.0098 


0.0340 


0.0340 


Qtot (b) 


-0.8627 


-0.8624 


0.4383 


0.4383 


0.3301 


0.3351 


1.1409 


1.1409 


ESin (MeV) 


1360.4385 


1360.4384 


1360.9970 


1360.9970 


1358.8995 


1358.8980 


1364.5696 


1364.5680 


Elin (MeV) 


827.3176 


827.3177 


827.3582 


827.3582 


827.1892 


827.1890 


828.1830 


828.1830 


E%t (MeV) 


-50.4839 


-50.4839 


-50.8174 


-50.8174 


-49.5856 


-49.5844 


-54.6531 


-54.6529 


E diT (MeV) 


365.7437 


365.7437 


365.7638 


365.7638 


365.7387 


365.7386 


365.9980 


365.9979 


Etot (MeV) 


-1024.7074 


-1024.7073 


-1024.6538 


-1024.6538 


-1024.3856 


-1024.3855 


-1023.4465 


-1023.4465 


J,, (h) 


0.00 


0.50 


0.00 


1.50 


0.00 


2.50 


0.00 


3.50 


J± {h) 


0.00 


0.82 


0.00 


0.00 


0.00 


0.43 


0.00 


0.01 



figurations and such a trial-and-error scheme, however 
helpful, is simply impossible to implement. Instead, we 
have resorted to a simple trick; namely, since the initial 
conditions do matter, we artificially generate slightly dif- 
ferent initial conditions for the two signature partners 
that converge to practically the same result. 

The main idea here consists in adding a tiny ro- 
tational frequency of about fta>=0.001MeV to break 
the degeneracy of signature configurations. This im- 
proves the convergence rate at the price of an in- 
significant numerical error of about 1-2 keV, on av- 
erage. We illustrate this fact by two specific exam- 
ples of blocked states in 163 Tb (native version, non- 
collective orientation) with different alignments. For 
the blocked state [411]l/2, at hu>=0 the total energy 
equals -1322.279268MeV, while for 7iu;=+0.001 and 
-0.001 MeV the total energies read -1322.279480 MeV 
(J|| = -n = -1/2) and -1322.279188 MeV (Jy = 
+Cl = +1/2), respectively. Similarly, for the blocked 
state [404] 7/2, the three corresponding energies are: 
-1321.725322 MeV, -132 1.725538 MeV (Jy = -ft = 
-7/2), and -1321.726010MeV (Jy = +Q = +7/2). To 
make our point, we deliberately show these energies with 
far too many digits than it is physically relevant. In these 
particular two examples, the numerical precision of the 
calculation is 10 eV, which is exceptionally good for odd 
nuclei. Therefore, the noted differences can only be at- 
tributed to the effect of the cranking term. Without this 
term, by blocking states /iq and jig one always obtains 
exactly the same HFB energy, E^ = Ep_ , because the 
full Skyrme functional is time-even. 



IV. RESULTS 

This section presents a number of Skyrme-HFB results 
for odd-mass nuclei. We begin by giving a detailed nu- 
merical comparison of the EFA with the exact blocking 
prescription in the limit of conserved time-reversal sym- 
metry. The impact of time-odd fields on the quasiparticle 
spectrum in the rare-earth region is shown in Sec. IIV Bl 
with the native, gauge, and Landau versions of SIII, SkP, 
and SLy4 functionals. The role of the nuclear alignment 
vector on physical observables is also studied. Results 
of calculations are compared with selected experimental 
data in Sec. IIV CI The triaxial polarization induced by 
a blocked quasiparticle is discussed in Sec. IIV Dl Finally, 
Sec. IIV El mentions the problem of the intrinsic insta- 
bility of certain Skyrme functionals that appear when 
time-odd terms are included. 



A. Validation of the EFA Approximation 

To demonstrate the numerical precision of our calcula- 
tions, TableUshows the results for four one-quasineutron 
states in 121 Sn obtained with hfbtho (efa) and HFODD 
(exact blocking). They are selected based on the mean- 
field spectrum ([To]) of 120 Sn. For the sake of this com- 
parison, the time-odd fields in hfodd have been switched 
off, thereby enforcing the regime where the exact block- 
ing procedure is strictly equivalent to the EFA, see Sec. 
Ill Dl Indeed, the obtained numerical differences between 
the efa and exact blocking are extremely small, less 
than 1 keV for the four different cases shown in Table 
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UJ This can be entirely attributed to various implemen- 
tations adopted differently in the two codes such as the 
method of computing the Coulomb potential, etc. The 
even-even core 120 Sn is spherical in its ground state. The 
quasiparticle blocking slightly polarizes the nuclear shape 
inducing small quadrupole deformations for some config- 
urations. 

Although the time-even observables obtained within 
the efa and exact-blocking are strictly identical if the 
time-odd fields are disregarded, this is not true for time- 
odd observables. In Table U this is illustrated by the 
values of alignments of the blocked quasiparticles aligned 
parallel (J||) or perpendicular (J±) to the symmetry axis. 
Of course, without time-odd fields, the direction of align- 
ment does not influence the time-even observables. 



the total energy is less than 1 keV for all of the excited 
states, and less than 0.002 b for the quadrupole moments, 
regardless of the quadrupole polarization exerted by a 
blocked state. 



B. Effect of Time-Odd Fields 

This section illustrates the effect of the various pre- 
scriptions for the time-odd channel © of the function- 
al. Calculations were performed for all odd-proton nu- 
clei with 63<Z<75 and 78<7V<104. For each of them, 14 
non-degenerate blocked configurations around the Fermi 
level have been considered. Altogether, 3,822 indepen- 
dent one-quasiproton states were studied. 



TABLE II: Comparison of EFA (hfbtho) with exact block- 
ing (hfodd) for ten one-quasiproton configurations in 163 Tb. 
Shown are the total energy Etot and total quadrupole mo- 
ment Qtot- The orbits are labeled by the asymptotic Nilsson 
quantum numbers [Nn z A]Q n . The SIII Skyrme functional is 
used in the particle-hole channel. Time-odd fields are disre- 
garded. The differences between hfbtho and hfodd results 
are shown in boldface. 





EFA ( 


'hfbtho) 


Exact 


(hfodd) 


Blocked State 


Qtot (b) 


1 Etot (MeV) 


Qtot (b) 


Etot (MeV) 


[411]3/2+ 


18.514 


-1323.495 


18.514 


-1323.495 


[532]5/2" 


17.759 


-1322.648 


17.759 


-1322.647 


[523]7/2~ 


18.554 


-1322.415 


18.555 


-1322.414 


[411] 1/2+ 


18.384 


-1322.322 


18.384 


-1322.322 


[413] 5/2+ 


18.654 


-1322.151 


18.654 


-1322.151 


[541] 1/2" 


20.138 


-1321.771 


20.136 


-1321.773 


[541]3/2" 


17.291 


-1321.357 


17.290 


-1321.357 


[530] 1/2" 


17.034 


-1320.762 


17.032 


-1320.762 


[420] 1/2+ 


17.766 


-1320.090 


17.767 


-1320.090 


[404] 9/2+ 


19.266 


-1319.851 


19.266 


-1319.851 



1. Native Functionals 

Table IIIII displays results for one-quasiproton states in 
163 Tb in the time-even, native, gauge, and Landau vari- 
ants of calculations. The alignment and elongation axes 
coincide with the y-axis of the reference frame. The time- 
even energies are shown in the absolute scale. For other 
variants, shown are displacements with respect to the 
time-even case: 



^^Todd _ ^Todd=0 



^Todd^O 



(31) 



In the particular example shown in Table IIIII the maxi- 
mum shift in one-quasiparticle levels due to the time-odd 
terms of the native functional is 127 keV. This is consis- 
tent with the earlier results of Refs. [7(| [7l| and over- 
all smaller than in the relativistic mean field approach, 
where time-odd polarization corrections seem to range 
from only a few dozen of keV in deformed actinides up 
to half a MeV in light nuclei (H, S3- 



Rare-earth nuclei provide an excellent testing ground 
for studies of deformed Nilsson orbits. Many of those 
nuclei are well-deformed, near-axial rotors and the de- 
formed mean field theory is particularly suitable to de- 
scribe their structural properties. Table ILTl shows a com- 
parison for several one-quasiproton configurations in a 
well deformed odd-proton nucleus 163 Tb. In hfbtho, 
the determination of a blocking candidate was made us- 
ing the mean-field spectrum (|16[) of the even-even core 
162 Dy. In the case of hfodd, to improve speed and sta- 
bility of the iterative process [49j], blocking candidates 
of a particle character (above the proton Fermi level of 
162 Dy) were selected from the mean-field spectrum of 
164 Dy while hole-like levels were selected from that of 
162 Dy. Of course, the final results do not depend on 
which particular even-even nucleus has been used as a 
core. 

The results shown in Table |TT] show again that without 
time-odd fields, the full blocking procedure is equivalent 
to the efa. It is worth noting that the difference on 



TABLE III: Energies (in MeV) of one-quasiproton config- 
urations in 163 Tb calculated using the time-even, native, 
gauge, and Landau variants of the SIII Skyrme functional, 
see Sec. Ill Bl Results for the time-even variant are shown 
in the absolute scale, while those for the other variants are 
shown as shifts (|31l) . 



Blocked State 


Time-even 


Native 


Gauge 


Landau 


[411]3/2+ 


-1323.495 


-0.075 


+0.042 


-0.125 


[532]5/2" 


-1322.647 


-0.052 


+0.029 


-0.105 


[523]7/2" 


-1322.410 


-0.060 


+0.039 


-0.080 


[411J1/2+ 


-1322.322 


-0.043 


+0.040 


-0.118 


[413] 5/2+ 


-1322.151 


-0.048 


+0.062 


-0.085 


[541] 1/2" 


-1321.773 


-0.007 


+0.055 


-0.075 


[541]3/2" 


-1321.357 


-0.047 


+0.036 


-0.107 


[530] 1/2" 


-1320.762 


-0.037 


+0.017 


-0.161 


[420] 1/2+ 


-1320.090 


-0.127 


+0.018 


-0.231 


[404] 9/2+ 


-1319.851 


-0.121 


+0.036 


-0.150 
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FIG. 1: Cumulative histogram of energy differences (|31|) for 
one-quasiproton states in the deformed rare-earth nuclei cal- 
culated with SIII, SkP, and SLy4 edfs. Dot-dashed open 
bins: native functionals; solid filled bins: Landau functionals. 
The bin size is 5keV. 



The overall impact of the time-odd fields on the energy 
of one-quasiproton states in the deformed rare-earth nu- 
clei is summarized in Fig. Q] which shows the distribution 
of AE Todd dSTJ) for SIII, SkP, and SLy4 edfs. When na- 
tive functionals are used, the total number of converged 
cases is 1,404 (524 for SIII, 443 for SkP, and 437 for 
SLy4). The average value of AE Todd is -50keV with a 
standard deviation of 42 keV. 

The magnitude of the odd-time effect depends on the 
choice of the edf. To illustrate this point, Fig. [5] dis- 
plays the distribution of Ai5 Todd for individual function- 
als. Focusing in this section on the native functionals 
(dot-dashed open bins), it is seen that the largest time- 
odd effect is predicted for SIII, which also shows an ap- 
preciable spread in values (configuration dependent). On 
the other hand, for the SkP parametrization the distribu- 
tion of AS Todd is fairly narrow, centered around -40 keV. 

By construction, Figs. [T] and [5] contain contributions 
from g.s. configurations and from nearly-lying excited 
states. Since g.s values are of particular importance as 
they impact mass predictions, Fig. [3] shows Ai? Todd for 
g.s. configurations only. The average value of the g.s. 
time-odd displacement is only ^50keV. Most of the few 
cases with |AE^ odd | > 150 keV correspond in fact to 
a collapse of pairing correlations in one of the 2 sets 
of calculations. It may be worth noting that the most 
recent HFB mass formula based on the Skyrme BSkl7 
parametrization yields a r.m.s deviation of 581 keV [lfjj . 
The uncertainty associated with neglecting the time-odd 
fields appears, therefore, to be smaller by an order of 
magnitude. 

In order to discuss the configuration dependence of the 
time-odd displacement, it is instructive to identify the 
s.p. orbits of interest. To this end, Fig. 0] shows the evo- 
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FIG. 2: Similar to Fig. [TJ except for individual Skyrme func- 
tionals: SIII (top), SkP (middle), and SLy4 (bottom). 



lution of the proton s.p. energies, defined as the eigen- 
values of the mean field operator (fl6|) . in the nucleus 
164 Dy calculated with SLy4 as a function of the axial 
quadrupole deformation 0:20 ■ This Nilsson diagram has 
been obtained by carrying out a set of constrained HFB 
calculations along a one-dimensional (Q20) path. 

Although the values of AE Todd are usually small, there 
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FIG. 3: Similar to Fig. \T\ except for g.s configurations only. 
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2. Landau junctionals 

Traditionally, only the time-even channel of Skyrme 
functionals has been adjusted to selected experimental 
data. That is, the time-odd channel has usually not been 
constrained. This is illustrated by the broad spread of the 
values of the isoscalar Landau parameters qg and g\ of 
the standard Skyrme functionals [U [72[ . In [23[ , a care- 
ful study of Gamow- Teller resonances within the Skyrme 
EDF theory yielded a set of 'optimal' Landau parameters 
that could be used to fix some of the coupling constants 
of the time-odd channel of the functional (namely the C t s 
and Cj). 

As seen in Tabic IIIIl the time-odd polarization in the 
Landau variant is greater than in the native variant, with 
the largest shift growing to 231 keV. The time-odd shifts 
in the gauge variant are generally smaller than for the 
native and Landau parameterizations. They also have 
opposite sign (time-odd polarization in the gauge variant 
decreases the binding energy while it is repulsive in native 
and Landau variants). 

The solid-filled bins in Figs. HH3] show A£: Todd for 
Landau-corrected functionals. The effect of this correc- 
tion is significant, as it shifts the centroid of most his- 
tograms by about 100 keV for SIII and SkP and 50keV 
for SLy4. When only ground-states are considered, the 
overall shift is of the order of 50 keV. 

To finish this section, let us recall that setting C^ s = 
was motivated in [23J to reproduce the energy and 
strength of the GT resonance, although different con- 
clusions about the role of this term were obtained later 
in [73|. In any case, the isoscalar channel governed by 
the Cq s term is not constrained by GT resonances, and 
in Refs. @, 0, [72| these terms are set to zero essen- 
tially to ensure the stability of the calculation. We briefly 
discuss this point in Sec. IIVEI 



FIG. 4: Nilsson proton levels of SLy4 in 164 Dy as function of 
the axial quadrupole deformation 0:20 («22 = 0). 



3. Alignments and Choice of the Quantization Axis 



are a few cases where the displacement can amount to 
more than 100 keV. A detailed analysis of the blocked 
configurations for all three interactions shows that the 
largest deviations correspond essentially to the [420] 1/2, 
[404]9/2, [400] 1/2, and [505] 11/2 Nilsson orbitals. It is 
interesting to note that the value of the s.p. angular mo- 
mentum does not seem to be crucial, since these states 
can be associated with both low and high-j spherical 
shells. In deformed rare-earth nuclei, equilibrium defor- 
mations are «20 ~ 0.25 — 0.30. As seen in Fig. @] in 
this deformation range, the orbital [420] 1/2 is a deep-hole 
state while [404]9/2, [400] 1/2, and [505] 11/2 are highly 
excited particle states. All these one-quasiproton exci- 
tations are strongly oblate-driving. A similar result has 
also been obtained for SLy4 and SkP. 



As discussed in Sec. Ill E[ one of the characteristic fea- 
tures of the treatment of odd nuclei in the blocking ap- 
proximation is the dependence of time-odd densities on 
the orientation of the alignment vector with respect to 
the principal axes or, equivalently, the choice of the self- 
consistent symmetries and quantization axis. To mea- 
sure this effect, we performed two sets of calculations. 
The first variant (_L) corresponds to the alignment vec- 
tor aligned along the y-axis, and the shape symmetry 
axis aligned along the z-axis. In the terminology of the 
cranking model, this case represents "collective rotation" 
perpendicular to the symmetry axis. In the second vari- 
ant (||), the nucleus is rotated by 90°, as described in 
Sec. MI Al so that the alignment and symmetry axes co- 
incide with the y-axis ( "non-collective rotation" ) . 

Note that in both situations y-signature and parity 
are conserved: the identification of blocking configura- 
tions via the position of the blocked state in a given sig- 
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FIG. 5: Cumulative histogram of deviations AE* = Eq P — Eq V 
between energies of one-quasiproton states calculated in the 
non-collective-rotation (||) and collective-rotation (_L) vari- 
ants. Dot-dashed open and solid filled bins correspond, re- 
spectively, to native and Landau variants. The bin size is 
5keV. 



nature/parity block hence provides a very robust way 
of tracking configurations before and after the Euler- 
rotation, as it is independent of the changes in other spa- 
tial characteristics of the quasi-particle wave-functions. 
As mentioned in Sec. IIII A[ the original Nilsson label of 
a q.p. state can be easily recovered in the non-collective 
orientation by simply exchanging the roles of the z and 
y axis in their computation. 

In Fig. [5] we show the distribution of differences AE^ = 

Eq P ~ E L for the 3,822 cases presented in the previous 
section (only those well converged are included in the 
plot). It is seen that the time-odd polarization due to 
the orientation of qp alignment gives an appreciable con- 
tribution to the time-odd shift, with the average value of 
AE^ being about 50keV in the native variant. The ori- 
entation effect seems to be weaker for Landau function- 
al. While the energy shift AE^ depends on the actual 
configuration, the total energy in the collective rotation 
scenario (_L) is overall lower than in the non-collective 
one (||) when native functionals are used. 



left panel), SIII (lower-left panel), and SLy4 (lower-right 
panel) functionals in the native variant. They are com- 
pared to experimental data. We follow the convention of 
Refs. [3(1 Lzij]j whereby the hole-like excitations are plot- 
ted below zero (representing the g.s. configuration) while 
the particle-like states are plotted above zero. 

The comparison with experiment suggests that the 
functional parametrizations employed in our work are 
not of spectroscopic quality for deformed nuclei. While 
the general deformation trends are reproduced and most 
of the orbitals found experimentally are indeed pre- 
dicted to appear around the Fermi level, the quantita- 
tive agreement with the data is not particularly impres- 
sive. For example, the SLy4 parametrization fails to re- 
produce the observed [523]7/2 g.s. of Ho isotopes; this 
state is predicted to lie 300-500 keV above the calculated 
[411] 1/2 ground state. Surprisingly, the oldest Skyrme 
parametrization SIII gives the best reproduction of ex- 
perimental band heads. The result of Fig. |5]is consistent 
with the conclusions of Ref. (38[; they found that the 
agreement of both spin and parity in the self-consistent 
models reaches about 40% for well-deformed nuclei re- 
gardless of the Skyrme force used. 

The three functionals used here have different isoscalar 
effective masses, m*=l, 0.707, and 0.7 for SkP, SIII, and 
SLy4, respectively. The effect of in* on shell structure 
is complex Q; among others, it impacts the density of 
states around the Fermi level. As seen in Fig. [6l the 
average level density obtained with SkP is indeed close 
to the experimental one. However, this does not nec- 
essarily mean that the spectroscopic properties are bet- 
ter described with this interaction: just as for SLy4, the 
ground-state is incorrectly assigned to the [411] 1/2 or- 
bital for all isotopes. 

There are, indeed, many factors that may impact the 
order of one-quasiparticle states. The recent analysis of 
spherical s.p. shell structure 0] has demonstrated that 
the isoscalar coupling constants in edf have a large im- 
pact on the position of s.p. energies and spin-orbit split- 
ting. It was also shown that the role of the effective-mass 
coupling constant cannot be reduced to merely changing 
the overall density of states. In fact, effective mass sig- 
nificantly influences relative positions of single-particle 
levels, including the splitting of spin-orbit partners. 



C. Experimental Odd-proton Spectra 

In well-deformed nuclei, one quasiparticle states can be 
related to the rotational band- head configurations [74j ]. 
In rare-earth nuclei, rich systematics of experimental 
data exist, and most importantly, the customary assign- 
ments of Nilsson labels [Nn z A]Q are available (30L ua]. 
Although these labels are approximate, they facilitate 
the comparison between theory and experiment. 

In Fig. |5] we show the one-quasiproton spectra for the 
Ho (Z=67) isotopic chain predicted with SkP (upper- 



D. Triaxial Shape Polarization 

Triaxial deformations of nuclear shape are enhanced 
at high spins [77], Lzl] • One spectacular example is the 
nuclear wobbling motion, which is caused by the fast ro- 
tation of triaxially-deformed nuclei [74], [79^ 81]. The phe- 
nomena of nuclear chirality is also tightly related to ax- 
ial asymmetry fH, H|[ . Recently, a systematic study 
of ground-state nuclear shapes in the framework of the 
macroscopic-microscopic model has also pointed to re- 
gions of triaxial instability in the nuclear chart |84j |. 

In the deformed rare-earth region that we consider in 
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FIG. 6: One-quasiproton band-head energies: calculated with SkP (top left), SIII (top right), and SLy4 (bottom left) Skyrme 
functionals, and extracted from experimental data (bottom right), plotted versus neutron number for Ho isotopes. Hole states 
are plotted below the ground state (zero energy) and particle states are plotted above. The time-odd terms are included in the 
native variant. Experimental data are taken from Refs. [3(il. Ft5(] . 



this work, the blocking of a quasiparticle built on intruder 
configurations has a strong 7-driving effect 85 88] . Most 
of the studies of this phenomenon are so far confined to 
high-spin states. Our calculations offer the opportunity 
to assess the degree of triaxiality in the g.s. configura- 
tions associated with weakly spin-polarized states. 

The calculated equilibrium deformations of one quasi- 
proton configurations considered in our survey are dis- 
played in Fig. [7J Time-odd terms are set to zero, so that 
results can be compared with the time-even calculations 
performed with hfbtho that define the axial reference 
point. As apparent in Fig. for the majority of con- 
figurations, triaxiality is very small, with 7 deformation 
typically less than 1°. 

Only a few highly excited states are characterized by 
a sizeable triaxial polarization: One such example is the 
state [402] 3/2, which originates from the spherical d 3 / 2 
orbital from the TV = 4 major shell and is pushed up into 



the TV = 5 major shell because of deformation. In Ta- 
ble IIVI we show the equilibrium deformations calculated 
with SLy4 for this specific configuration in a number of 
isotopes. The excitation energies of [402] 3/2 range from 
0.2 to 3MeV. On average, the net energy gain induced 
by the triaxial polarization of the core is of the order of 
200 keV in this extreme case. 

As indicated, the results presented in Fig. [7] have been 
obtained by setting all time-odd fields to zero. When this 
constraint is released, g.s. configurations remain over- 
whelmingly axial, independently of the orientation of 
the alignment vector, cf. discussion in Sec. Ill El How- 
ever, we do observe that in the collective orientation 
limit, low-j intruder states such as [541] 1/2 (from h.9/2) 
and [550]l/2 (from h n / 2 ), or high-j intruder states such 
as [505]ll/2 (from hn/2), seem slightly more unstable 
against 7-polarization than in the non-collective situa- 
tion. 
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FIG. 7: Distribution of the equilibrium deformations in the 
(X,Y) plane, X = /3cos( 7 + tt/6) and Y = /3sin( 7 + tt/6), 
where /3 and 7 are the standard Bohr quadrupole deforma- 
tion parameters. The figure corresponds to the 3,822 differ- 
ent blocked configurations considered in this study. Time-odd 
terms are set to zero. 



TABLE IV: Equilibrium deformation of the [402]3/2 blocked 
configuration in several odd-proton isotopes with the SLy4 
interaction (time-odd terms set to zero). A 1 E represents the 
gain (in keV) in energy induced by the triaxiality. 



z 


N 


E* (MeV) 


P 


7 (deg) 


A 7 £ (keV) 


69 


90 


1.506 


0.21 


7.7 


-191 


69 


92 


2.070 


0.25 


6.7 


-187 


69 


94 


2.471 


0.28 


5.9 


-191 


69 


96 


2.745 


0.29 


5.4 


-184 


69 


98 


2.955 


0.30 


4.0 


-124 


71 


86 


0.232 


0.13 


19.6 


-233 


71 


88 


0.647 


0.17 


11.8 


-214 


71 


90 


1.106 


0.20 


8.9 


-195 


7:i 


88 


0.442 


0.16 


8.1 


-203 


73 


90 


0.717 


0.18 


8.9 


-205 



This overall axial stability is illustrated in Fig. HI where 
the distributions of the 7 angles for well-deformed odd- 
proton states in the rare-earth nuclei are plotted. For 
better legibility of the figure, the very rare pronounced 
triaxial cases with 7 > 2° have been omitted - they have 
been discussed above, and so have the many near-axial 
states with 7 < 0.05°. In the lower panel corresponding 
to the collective orientation, the few points beyond 7 = 
1° correspond to the 7-driving orbitals. If the rotational 
frequency ui y is increased (cranking), we find that the 
degree of triaxiality increases accordingly [H, H3, HH ■ 
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FIG. 8: Triaxiality (measured by the angle 7) of well- 
deformed odd-proton states (/3 > 0.1) in the rare-earth region 
calculated with SIII, SkP, and SLy4 edfs for the two orien- 
tations discussed in Sec. IIVB 31 collective (upper-panel) and 
non-collective (lower panel). 



E. Finite-size Instabilities of Band-head 
Calculations 



It has been shown that some parametrizations of 
the Skyrme ener gy functional could be prone to finite 
size instabilities [90-92] . In particular, the time-even 
and time-odd C As s f As t terms could, in some 



cases, lead to divergences of the HFB iterative procedure. 
The size of these instabilities depends on a number of 
factors such as the edf parametrization, particle num- 
ber, and specific implementation of the DFT solver. The 



detailed analysis of edf instabilities performed in [92 1 
has been based on the RPA response function approach 
of Ref. [93| implemented to Skyrme functionals |94l [95j j . 
Results were reported in 40 Ca and 56 Ni for the SkP and 
SLy5 parametrizations. 

Finite-size instabilities governed by C t As terms are am- 
plified in polarized systems such as odd-mass nuclei. In- 
deed, these terms are only active when time-reversal sym- 
metry is broken. As was shown in Sec. lIVBl the impact of 
time-odd components is weak, at least in the rare-earth 
region that we study. It is therefore possible to scale 
these terms by slightly varying the values of C As , with- 
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out impacting significantly the calculated properties. By 
contrast, scaling the coupling constants C t Ap could result 
in totally non-physical solutions. 

According to 92], the functionals employed in this 
work, namely SIII, SkP, and SLy4, should not be par- 
ticularly sensitive to spin instabilities. Indeed, the rate 
of convergence in our calculations is of the order of 40- 
50% for those three cases. This is less than for even-even 
axially deformed nuclei, but such a rate can be tied to 
factors such as collapse of pairing, level crossings, etc. 



70 - 
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Cf /q s (std) 

FIG. 9: Convergence rate of hfb equations with SkP, SkO, 
and SkM* functionals for one-quasiproton states in odd- A Ho 
isotopes with 88<iV<104 as a function of the scalar-isoscalar 
coupling constant Cq. See text for details. 

However, other Skyrme parametrizations may be prone 
to severe and systematic divergences. To illustrate this 
point, we have performed a set of calculations with three 
functionals: SkO 0, SkP, and SkM* [97]. For each 
of those, we have used the native variant of the time- 
odd terms; only Cq s is multiplied by a scaling factor a 
ranging between zero (no coupling) and one (standard 
coupling) . A measure of stability of the iterative process 
is the rate of convergence for a pre-defined set of one- 
quasiparticle states. A result is deemed converged if the 
binding energy does not change by more than 2 keV from 
one iteration to the next for 3 consecutive iterations. We 
show in Fig. [S] the evolution of this convergence rate as 
a function of a. Our set of configurations consists of 24 
different one-quasiproton states in nine odd-A Ho iso- 
topes with 88<A^<104. Therefore, the sample size used 
to define the convergence rate is 216. 

According to Fig. [9j SkM* and SkP parametrizations 
are stable with respect to variations of Cq, but the SkO 
functional exhibits a sharp drop in the convergence rate 
when a > 0.5, i.e., C^ s > 35MeV. Preliminary inves- 
tigation of the RPA response function [§8| suggests that 
instabilities could occur for transferred momenta q of the 
order of 2.2-2.5fm~ 1 for this particular value of C^ s . 



These results nicely agree with the original findings of [92] 
and emphasize the need of testing edfs against finite-size 
instabilities. 



V. CONCLUSIONS 

In this work, we carried out the systematic theoret- 
ical survey of one-quasiproton states in deformed rare- 
earth nuclei. Our study is based on the symmetry- 
unconstrained Skyrme hfb framework that fully takes 
into account time-odd polarization effects. 

We show that the equal filling approximation is equiv- 
alent to the full blocking when the time-odd fields are 
put to zero. In this case, an arbitrary combination of 
time-reversed orbits can be used to define the blocked 
orbit, and this can be nicely quantified by introducing 
the notion of alispin. We emphasize the role of symme- 
tries, and in particular nuclear alignment properties, in 
the exact treatment of the blocked state. 

Our systematic survey indicates that, when native 
functionals are employed, the contributions from time- 
odd fields to the energy of the g.s. and low-lying excited 
states is rather small, around 50keV on average, with a 
variation around 100-150 keV. Significant differences are 
found from one interaction to another, although the ef- 
fect remains small for the three interactions considered. 
Correcting the time-odd channel (Landau functionals) in- 
creases the contribution of the time-odd channel to the 
total energy by about 50%. For the functionals in the 
gauge variant, the time-odd effects are weak and oppo- 
site in sign. 

By explicit calculations we demonstrated that the 
choice of the alignment orientation with respect to the 
quantization axis does impact predicted time-odd polar- 
ization energies. The resulting energy shifts are appre- 
ciable in the scale of predicted time-odd displacements. 

Standard paramctcrizations of the Skyrme interaction, 
such as the SIII, SkP, and SLy4, give a qualitative, 
but not quantitative description of experimental one- 
quasiproton spectra in the rare-earth region. We find 
that the triaxial shape polarization effects are generally 
small in the nuclei considered. Finally, we point to the 
sensitivity of DFT calculations for one-quasiparticle states 
to finite-size instabilities of the underlying edf. A de- 
tailed investigation of this effect is currently under way. 

The weak impact of the time-odd fields on spec- 
troscopic properties implies that global studies with 
symmetry-restricted HFB solvers such as hfbtho could 
be very useful to extract information related to the 
isovector properties, shell structure, and shapes. For such 
a purpose, time-odd fields may be safely neglected. 
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